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Abstract. We simulate a cluster of 144179 stars including 13107 primordial 
hard binaries (10% of the total number of single stars and binary centers of mass), 
using direct integration of the equations of motion of all stars and binaries and 
incorporating the effects of stellar and binary evolution. The initial conditions 
are representative of young dense star clusters in the Local Group and other 
nearby galaxies like the Antennae and M82. We find that the early phase of 
core collapse, driven by mass segregation, is not appreciably delayed by the 
presence of a large number of hard binaries. By the end of the simulation, at 
an age of 115 Myr, the cluster radius has expanded by about a factor of two. 
This may be explained as adiabatic expansion driven by the loss (via stellar 
evolution) of ~ 40% of the initial total mass. Binary dynamics apparently has 
little effect on the early cluster expansion. During the evolution, the total binary 
fraction drops at a roughly constant rate of ~ 0.01% per Myr. The fraction of 
very hard binaries, however increases at about 0.025% per Myr. By the end of 
the simulation the cluster contains 37 binaries containing at least one black hole; 
roughly half (17) of these contain two black holes. 



Introduction 



Young and dense star clusters (YoDeCs) are a relatively new class of object, 
containing about 10 4 -10 6 stars and having ages less than about 10 Myr and 
densities exceeding a few times 10 5 stars/pc 3 . There are four examples in the 
Milky Way Galaxy: the Arches and Quintuplet clusters near the Galactic center, 
and NGC 3603 and Westerlund 1 in the disk. Tabled summarizes the observed 
parameters for some known YoDeCs. In addition to the clusters just mentioned, 
we also include R 136, the central cluster in the 30 Doradus region, and MGG 11, 
in the starburst galaxy M82, which may contain an intermediate-mass black hole 
( Portegies Zwart et al.l|2004D . 

Motivated by these observations, as well as by theoretical considerations, 
we distinguish two families of YoDeC: those which are effectively isolated, and 
those which are strongly perturbed by the tidal field of their parent Galaxy. 
We currently know only two clusters in the latter category — the Arches and 
Quintuplet systems. The other clusters in Tabled are relatively unperturbed. 
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In this paper we focus on the dynamics of the primordial binary population in 
unperturbed young dense clusters. 

The binary frequency and the distribution functions of binary orbital pa- 
rameters are poorly constrained by observations. For two YoDeCs we have a 
rough estimate of th e binary fraction: NGC 3603, with an estimated binary 
fraction of about 0.3 ( Stolte et al.l l200"4T ) . a nd R136, which may conta in up to 
70% binaries among its most massive stars ( Portegies Zwart et al.ll200'2T ) . These 
indications of binary frequencies are based on the X-ray signatures of these clus- 
ters. Further evidence for a hig h binary fraction in two star clusters near the 
Galactic center is presented by ( Yusef-Zadehll2003T ) . 



2. The simulation model and initial conditions 

As initial conditions for our simulations we use a cluster of 131072 stars drawn 
from a Salpeter initial mass function between 1 and 100 M Q . The lower mass 
limit is motivate d by the observed age and mass-to-light ratio of the young star 
cluster MGG11 riMcCradv et al.ll2003D . A randomly selected subset of 13107 



(10%) stars receives a secondary companion with mass between IMq and the 
mass of the selected (primary) star. Orbital separations range from Roche-lobe 
contact (at about E ~ 10 3 kT, where |fcT is the mean stellar kinetic energy) to a 
binding energy of E = 10 kT (at an orbital period of about a week), while initial 
eccentricities are assumed to follow a thermal distribution (see the left panel 
in Figure0J. The initial half-mass radius of our simulation is 3.2 pc, and even 



Table 1. Observed properties of selected young, dense star clusters, mostly in 
or near the Milky Way. Mass (M), age and half-mass radius (thm) are taken from 
the literature. The last two columns give the suspected binary fractions and an 
estimate of the initial mass function. Two question marks for the binary fraction 
indicates that binaries are likely to be present but it is impossible at this stage 
to give a reliable estimate, but within the next few years better estimates may 
be expected. A horizontal dash indicates that nothing is known about the binary 
fraction in that cluster. The mass function is designated as a Power-law (PI) with the 
exponent between parenthesis, or as a Salpeter function (identified by S), a power- 
law with exponent 2.35, with lower and upper limits in Mq as sub- and superscript, 
respectively. 
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References: a) Figer et al. (1999;2002); b) Glass et al. (1987); c) Brandl (1999), ; d) Vrba et 
al. (2000); e) Brandl et al. (1996), ; f) McCrady et al. (2003) 
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though our simulated cluster is isolated we choose a King model with Wo = 12 
as the initial density profile. 

We solve the equations of motion of all the stars in the cluster, at the 
same time calculating the evolution of the stars and b inaries using the Star- 
lab 1 software environment (iPortegies Zwart et al.ll200lh . We use the GRAPE- 
6 special-purpose computer (|Makino et al.l 120031 ) to speed up the calculations. 
Even though no older siblings of the YoDeCs listed in Tabled are known, our 
cluster simulations generally continue to ages of about 100 Myr. 

Rather than show the usual evolution of the core and Lagrangian radii we 
show instead in Figure^^c simulated images of the model cluster, one near the 
start of the simulation, and one nea r the end. Clearly, the cluster expands con- 



siderably with time. As discussed in ()Portegies Zwartl2 004h the global evolution 



is characterized by three quite distinct stages: A an early relaxation-dominated 
phase, followed by phase B, during which the ~ 1 % (by number) most mas- 
sive stars evolve quickly and lose an appreciable fraction of their mass. Finally, 
in phase C, stellar evolution sl ows and relaxation takes over until the cluster 
dissolves (|Portegies Zwartll2004 ). For this particular simulation, phase A lasts 
for about 3 Myr, phase B lasts until ~ 50 Myr, followed by phase C for the 
remainder of the run. 
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Figure 1. Simulated HST NICMOS images of our model cluster at a dis- 
tance of 8kpc, at birth (left) and at an age of about 115 Myr (right). Images 
are taken with identical observation time in the J-band; the image field of 
view is about 10 pc across. 



3. Results 

Figure|2] shows the evolution of the number of binaries in our simulated cluster. 
The total number of binaries drops by ~ 9% in 115 Myr, whereas the number of 



x Starlab is available at http://www.manybody.org/starlab/starlab.html 
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stars drops by only 1%. The total cluster mass drops by about 40% during the 
same period. 

The top solid curve gives the total number of binaries throughout the simu- 
lations, from birth to an age of 115 Myr. The lower solid curve gives the number 
of binaries within 10 r c . We use ten times the core radius (r c ) here because the 
number of binaries withing one core radius is so small. The dashed curve shows 
the number of very hard binaries, having binding energy E > 100 kT. 




50 100 



t [Myr] 

Figure 2. Evolution of the number of hard binaries including binaries softer 
than 10 kT for our simulation (top solid curve). Binaries with binding energy 
> 100 kT are presented as the dashed curve, and the number of hard bina- 
ries within 10 r c are indicated by the middle solid curve. The lower solid 
curve shows the binaries containing at least one compact object — black hole, 
neutron star or white dwarf. An enlargement of this curve is presented in 
Fig.El 



In Figure|3]we present, as functions of time, the number of binaries with at 
least one compact object, black hole, neutron star or white dwarf. The number 
of black holes in binaries rises sharply shortly after the start of the simulation 
with a peak near 8 Myr. This is the moment when the turn off drops below 
the minimum mass (~ 20 Mq) for forming black holes: lower-mass stars form 
neutron stars. This threshold is also visible as the sharp increase of the number 
of binaries containing a neutron star. The number of black holes in binaries 
drops rapidly from this moment on, because many of their companions form 
neutron stars in supernova explosions. These newly formed neu tron stars receive 
a mu ch higher asymmet ric kick velocity during t heir formation ( Lvne Lorimerl 
I1994T ) than black holes dGualandris et al.1 12004 ). Note that binaries containing 
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a neutron star and a black hole are counted twice in this figure, once among the 
(bh, *) and once for (ns, *). 

The middle solid curve in Figure|3] shows an interesting dip between about 
75 Myr and 80 Myr binaries within 10 r c are depleted. The explanation for this 
local depression is the occurrence of core collapse. This is rather surprising, as 
the energy stored in primordial binaries is sufficient to support the cluster core 
and to prevent core collapse. 

For a star cluster with a broad mass spectrum and without primordial 
binari es, core collapse occurs in about 2 0% of the initial half-mass relaxation 
time (Mc Millan &: Portegies Zwartl 120 03). With a half- mass relaxation time of 
~ 360 Myr in our simulation, one therefore expects core collapse in about 70 Myr. 
We attribute the fact that the observed collapse happens about 10 Myr later than 
that to the effect of stellar mass loss, as by this time the cluster has lost about 
36% of its mass by stellar evolution. (Note that less mass would have been 
lost by the cluster if we had adopted a lower limit to the initial-mass function.) 
Most of the mass loss happens in the cluster core, effectively driving an adiabatic 
expansion; it appears that binary heating does not play an important role in the 
core evolution. 

Even though we have not presented here definitive proof that the binary 
population is ineffective in influencing this phase of the core evolution, the col- 
lapse of the cluster core suggests that binary heating delays core collapse by at 
most 10 Myr. This may be contrasted with the results of earlier Monte-Carlo 
simulations, which indicate that core collapse is si gnificantly delayed by the 
presence of a rich population of primordial binaries ( Fregeau et al.ll2003h . 

The first white dwarfs appear at 23.8 Myr, at a turn-off mass of about 
10 M©. Isolated stars (i.e. unperturbed by any companion) are expected to evolve 
into white dwarfs if their initial mass is less than roughly 8M Q . However, in 
binary systems, early stripping of the hydrogen envelope may cause more massive 
stars also to become white dwarfs instead of collapsing to neutron stars. The 
population of compact binaries of clusters older than about 40 Myr is dominated 
by white dwarfs, although at a formation rate of about 100 white-dwarf binaries 
per 30 Myr, it will take about 2 Gyr before they become more common than 
primordial main-sequence binaries. 

Figure^] shows the distribution in semi- major axis and eccentricity of the 
binary population at the start of our simulation (left panel) and at 115 Myr 
(right). We have rendered the image in gray-shades, with darker shades indicat- 
ing more systems within the bin. The left panel in Figure0]represents essentially 
the zero-age population, except that we show the system after about one cluster 
crossing time (0.1 Myr) to eliminate (for example) any spurious initializations of 
eccentric binaries which should circularize shortly after formation. The orbital 
separation ranges from about 15 R0 to about 1000 Rq. The lack of very high 
eccentricities (e ^, 0.9) is mainly due to tidal effects within binary systems. An 
encounter is unlikely to ionize a binary, but effectively changes its eccentric- 
ity; highly eccentric binaries (e ^ 0.6) typically become more circular after an 
encounter. 

In the right panel of Figure0]we show the distribution in semi-major axis 
and eccentricity for the binaries remaining in our simulation at an age of 115 Myr. 
The two most striking differences compared to the left (zero age) panel are the 
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Figure 3. Magnification of the evolution of the number of compact objects 
in the simulation. The solid curve shows the number of binaries containing 
at least one black hole, neutron star or white dwarf (identical to the bottom 
solid curve in Figure|2J) . The dashed curve gives the number of binaries with at 
least one neutron star, indicated as (ns, *). The two dotted curves, indicated 
as (bh, ★) and (wd, *), give the number of binaries containing at least one 
black hole and at least one white dwarf, respectively. 
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Figure 4. Distribution of the orbital separation (in Rq) and eccentricity 
of all binaries in our simulation near zero age (left frame) and at the end 
of the calculation, at about 115 Myr (right frame). Darker shades indicate 
more binaries, as is indicated in the normalization in the upper-left corner of 
the left panel. The total number of binaries and the age of the cluster are 
indicated near the top of each frame. 
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depletion of binaries with orbital separations a ^ 200 Rq, and the enormous 
increase in binaries on circular orbits. 

The pile-up in Figure0]of wide (a ^ 200 Rq) binaries with circular orbits and 
the depletion of wide eccentric binaries is caused mainly by the internal evolution 
of the binaries in combination with the initial conditions. The initial binary 
separations were selected between Roche-lobe contact and a binding energy of 
E = 10 kT. The hard-soft boundary for a binary with orbital separation a and 
total mass m may be estimated as 



Here N, M, and R are the number of stars and the total cluster mass and 
half-mass radius. Binaries with total mass exceeding 5Mq will, according to 
our initial requirement that all binaries have binding energy E > 10 kT, have 
orbital separations smaller than 200 Rq . The turnoff mass for our cluster at an 
age of 115 Myr is about 4.5 Mq, so binary components more massive than that 
are either ascending the giant branch or have already done so. The binaries 
with initial orbital separations exceeding ~ 200 Rq host mainly rather massive 
stars, which have already evolved off the main sequence by the time the cluster 
is 115 Myr old, and therefore these binaries have experienced strong internal 
evolution. Note here that wide binaries tend to shrink upon Roche-lobe over- 
flow, due to angular momentum loss via non-conservative mass transfer, whereas 
shorter-period binaries tend to expand; this leads to a pile-up of circular bina- 
ries around a = 100-300 Rq. Binaries containing lower- mass stars generally have 
tighter orbits, and do not populate this region, at least not during the cause of 
the simulation. 

The few eccentric binaries with orbital separations a ^ 200 Rq at an age 
of 115 Myr have experienced two supernovae; in most cases (11 out of 14) both 
supernovae led to the formation of a black hole. The few binaries to the left of 
the circularization boundary are in the process of rapid circularization due to 
tidal effects. 

4. Summary 

We report the results of a large (128k) direct N-body simulation with 10% hard 
(E XL 10 kT) primordial binaries. Our simulations include the effects of stellar 
mass loss, collisions, and coalescence, and the internal evolution of binaries. 

Our main conclusions may be summarized as follows: core collapse in the 
simulation occurs at an age of about 80 Myr. In a star cluster with similar 
initial conditions but without stellar mass loss and without primordial binaries, 
core collapse would be expected to occur at about 70 Myr. We therefore argue 
that core collapse is slightly delayed by stellar mass loss and the presence of the 
binaries. The former appears to be the main cause of the delayed core collapse. 

The number of binaries containing one or two compact objects increases 
steadily from the start of the simulation. After about 30 Myr, the number of 
binaries with at least one black hole remains roughly constant at about 40, 
whereas the number of binaries containing at least one white dwarf continues to 
grow at a roughly constant rate of about 100 per 30 Myr. 
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